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Abstract 

This paper establishes limit theorems for a class of stochastic hybrid 
systems (continuous deterministic dynamic coupled with jump Markov 
processes) in the fluid limit (small jumps at high frequency), thus extend- 
ing known results for jump Markov processes. We prove a functional law 
of large numbers with exponential convergence speed, derive a diffusion 
approximation and establish a functional central limit theorem. We apply 
these results to neuron models with stochastic ion channels, as the number 
of channels goes to infinity, estimating the convergence to the determinis- 
tic model. In terms of neural coding, we apply our central limit theorems 
to estimate numerically impact of channel noise both on frequency and 
spike timing coding. 
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1 Introduction 

In this paper we consider stochastic hybrid systems where a continuous de- 
terministic dynamic is coupled with a jump Markov process. Such systems 
were introduced in [6] as piecewise deterministic Markov processes. They have 
been subsequently generalized to cover a wide range of applications: commu- 
nication networks, biochemistry and more recently DNA replication modeling 
[21 [T5l [20l [23] . We are interested in the fluid limit for these systems considering 
the case of small jumps of size 1/N at high frequency N, with a view towards 
application to neural modeling. 

The general class of model we consider is described in section 2.1, and for 
the sake of clarity, we describe here a simple example which retains the main 
features. Consider a population of N independent individuals, each of them 
being described by a jump Markov process Uk(t) for k = 1, ...,N with states 
and 1, and with identical transition rates a > 0,f3 > as follows: As an 




empirical measure, we define the proportion of individuals in state 1 at time t 
by: 

1 N 

ejv(*) = J^^2 u k(t) 
fe=i 

The model becomes hybrid when we assume a global coupling through a variable 
V/v € R, in the sense that the rates a(V/v) and /?(V/v) are functions of Vjv- This 
variable Vn is itself solution of a differential equation, between the jumps of 
e N (t): 

dVN tnr \ 

where / : R 2 — > R. In the general case, this model is extended with more 
general non-autonomous jump Markov processes, the global variable can be 
vector valued and the transition rates can be functions of the empirical measure 
(section 2.1). 

We prove convergence in probability on finite time intervals, with techniques 
inspired by p], of the solution of the stochastic hybrid system to a deter- 



2 



ministic limit x = (v,g). For the example above, x is solution of: 



dv 
It. 
dg 



f(v,g) 



(1 - g)a(v) - gp(v) 



We derive a diffusion approximation and prove a functional central limit theorem 
that helps characterizing the fluctuations of both the discrete and continuous 
variables around the deterministic solution. We obtain that these fluctuations 
are a gaussian process which corresponds to the asymptotic law of the linearized 
diffusion approximation. We further obtain an exponential speed of convergence 
which relates the tail distribution of the error En{T) — s\ip\X^ — x\ 2 to the 

[0,T] 

size parameter N and the time window T : for A > and N large, 



Thus the convergence result can be extended to large time intervals [0,T(N)], 
provided that T = T(N) is such that NH(T(N)) ->■ oo. Inequality (1.1) is a 
new result which provides an estimate to the required number N of individuals 
to reach a given level of precision. This number increases with the time scale on 
which one wants this precision to be achieved. For system subject to finite-size 
stochasticity, sometimes called demographic stochasticity it provides a relation 
between the reliability time-scale to the population size N. There are other 
ways of obtaining a law of large numbers, for example using the convergence of 
the master equation or of the generators [TU]. We want to highlight here that 
our proof is based on exponential inequalities for martingales. Other ways of 
obtaining a law of large numbers would not be likely to provide an estimate 
such as (1.1). 

Our mathematical reference on the fluid limit is the seminal paper 22 which 
contains a law of large numbers and a central limit theorem for sequences of 
jump Markov processes. Recently, a spatially extended version of these models 
has been considered in [T], for a standard neuron model. The author shows 
convergence in probability up to finite time windows to a deterministic fluid 
limit expressed in terms of a PDE coupled with ODEs. In the present paper, we 
consider a class of non-spatial models which however includes multi compart- 
mental models, by increasing the dimension. We extend the results of [22] to 
stochastic hybrid models at the fluid limit. 

Neurons are subject to various sources of fluctuations, intrinsic (membrane 
noise) and extrinsic (synaptic noise). Clarifying the impact of noise and vari- 
ability in the nervous system is an active field of research [29], [11]. The in- 
trinsic fluctuations in single neurons are mainly caused by ion channels, also 
called channel noise, whose impacts and putative functions are intensively in- 
vestigated (37] |30j [28] , mainly by numerical simulations. Our motivation is to 
study the intrinsic fluctuations in neuron models and we think that stochastic 
hybrid systems are a natural tool for this purpose. The channels open and close, 



P(E N (T) > A) < e 



ANH(T) 



(1) 
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through voltage induced electromagnetic conformational change, thus enabling 
ion transfer and action potential generation. Because of thermal noise, one of 
the main features of those channels is their stochastic behavior. 
In terms of modeling, our starting point is the stochastic interpretation of the 
Hodgkin- Huxley formalism [TB]. In this setting, ion channels are usually mod- 
eled with independent Markov jump processes, whose transition rates can be 
estimated experimentally [35]. These stochastic discrete models are coupled 
with a continuous dynamic for the membrane potential, leading to a piecewise- 
deterministic Markov process. Thus, the individuals are the ion channels and 
the global variable Vn the voltage potential (cf. section 3.). Deterministic hy- 
brid kinetic equations appear to be a common formalism suitable for each stage 
of nervous system modeling as shown in [8] . This latter study provides us with 
a framework to introduce stochastic hybrid processes to model action poten- 
tial generation and synaptic transmission, as stochastic version of deterministic 
kinetic models coupled with differential equations through the transition rates. 
On the side of neuron modeling applications, the limit behavior of a similar but 
less general model is considered in |12j , using an asymptotic development of the 
master equation as N — > oo, which formally leads to a deterministic limit and a 
Fokker-Planck equation (Langevin approximation) , providing the computation 
of the diffusion coefficients. The Langevin approximation is also studied in |34j . 
but in a simplified case where the transition rates are constants (independent of 
Vn), which is actually the case studied in [32]. Our mathematical results extend 
these previous studies to a wider class of models (if we put aside the spatial as- 
pects in p]), providing a rigorous approach for the Langevin approximation, and 
establishing a central limit theorem which describe the effect of channel noise 
on the membrane potential [32j . The convergence speed provides a quantitative 
insight into the following question : if a neuron needs to be reliable during a 
given time-scale, what would be a sufficient number of ion channels? We thus 
provide a mathematical foundation for the study of stochastic neuron models, 
and we apply our results to standard models, quantifying the effect of noise on 
neural coding. In particular, both frequency coding (sec. 3.5.1) and spike tim- 
ing coding (sec. 3.5.2) are numerically studied with Morris-Lecar neuron model 
with a large number of stochastic ion channels. 

Generically, stochastic hybrid models in the fluid limit would arise in multiscale 
systems with a large population of stochastic agents coupled, both top-down 
and bottom-up, through a global variable, leading to an emergent cooperative 
behavior. Starting from a microscopic description (ion channels), the central 
limit theorem as stated in this paper leads to a description of the fluctuations of 
the global variable (membrane potential). So, in the perspective of applications, 
it would be interesting to investigate how our framework and results could be 
developed in other fields than neural modeling: for instance in chemical kinetics, 
in population dynamics, in tumor modeling, in economics or in opinion dynamics 
theory. In a more mathematical perspective, it would be interesting to consider 
a wider class of models, for instance by including spatial aspects as in [I] or by 
weakening the independence assumption. Other questions could be investigated, 
for instance concerning escape problems, first passage times and large deviations, 
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whenever TV is large or not. 

Our paper is organized as follows. In section 2. we define our model and 
formulate the main results. In section 3., we apply our results to neuron models. 
In section 4. we give the proof of the law of large numbers and its convergence 
speed (Theorem 2.1) and in section 5. we give the proof of the Langevin 
approximation (Theorem 2.2) and central limit theorems (Theorem 2.3- 
2.4-2.5). 

2 Model and main results 

This section contains the definition of our general model and states the main 
theorems. 

2.1 Model 

Stochastic hybrid model (S N ) Let p,q,N e N*, and rj e N* for all 1 < 

j < q. Let d — Y^j=i r j- We define the stochastic hybrid model (SW), whose 
solution 

X N (t) = (V N (t),e N (t)) e R p x R d , t > 

satisfies: 

and e 7v = (e^\...,e^) with e R rj , where the processes e${t) are q in- 
dependent jump Markov processes. Note that the differential equation for Vn 
is holding only between the jump times of the process ejv, with updated initial 
conditions. For 1 < j < q, processes e^\t) are characterized by, 

• their state space : E$ = {(xi, x Tj ) € {0, jf, 1}^ \ YJk=i x k = ^} 

• their intensity A$: for X = (V, e) e R p x R d , \$(X) = NX^(X) with 

fe=i i=i, ijtk 



their jump law /j,^: we define «o' = (0, 0, 1, 0, 0) G R rj and 



b 



(i) (i) 

Ua — ul for 1 < a, b < rj. The transition of an individual agent in the 
population j from one state a to another state b corresponds to a jump of 
z = jjU^ a for the process e$ . Thus we define: 



X + -AXi tb = (y,eW,...,eV),..., e M) 



gO) = e (i) + u « 
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So that the jump law for a jump of z is given by: 

U) (x z) _ Ca^aA*) if = i 
Mjvl ' j " AO)(X) N b > a 



for all 1 < a, & < r*j such that ea 7^ and el ^ 1, and 
fJ^(X, z) — otherwise. 
For a more formal definition we refer to [5J. 

For 1 < k < rj, the fc-th component {e$}k of vector ejy can be interpreted as 
the proportion of agents of type j which are in the state k in a population of 
size N. 

We show below in Theorem 2.1 that this stochastic hybrid model has a limit 
as N — ► 00 which is the following deterministic model. 



Deterministic model (D) We define the deterministic model (D), whose 
solution X = (v,g) e R p x R d with g = (gW, ■•,g (9 - ) ) satisfies: 

( v = f(v,g) 

l l<j:<rj,i^fc 

for all 1 < j < g, VI < fc < r 3 . The first equation is the same as in the stochas- 
tic model (deterministic part) and the second equation corresponds to the usual 
rate equation, with a gain term and a loss term. 

The following example illustrates the general model in a simpler relevant setting 
motivated by applications. This setting will be used in the proofs in section 4 
and 5 in order to make the arguments clearer. 



Example We consider the case where p = q = 1 and n = 2. We can construct 
a stochastic hybrid process as follows: first let us introduce a collection of N 
independent jump Markov processes for 1 < k < N with : — > 1 with 
rate a(Vjv) and f — > with rate /3(V/v) : 

<>(Vn> 

p-i 1 



where Vjv is defined below. We then consider ejv(t) = ({ejv}o(*), { e Jv}i(i)) the 
proportions of processes in the states and 1 . In this case, the stochastic hybrid 
model (Sn) can be written as: 



G 



' V N (t) = f{V N {t),e N {t)) 
< e N (t) = (±Eti6o(4%±EL^M k) )) 
, V N (0) = v : ejv(O) = (u , 1 - u ) 

Note that, if we define UN{t) — j; J2k=i ^i( u t^)i tnen we have ejv(i) = (1 — u N(t), u N 
so that the solution is determined by the pair Xff(t) = (V/v(t), UN(t)). 
Thus, each member of the sequence of jump Markov processes {un}n>i is 
characterized by 

• its state space E N = {0, 1}, 

• its intensity \ N (V N (t),u) = N[u/3(V N (t)) + (1 - u)a(V N (t))]. This inten- 
sity is time-dependent through V/v(t). 

• its jump law 

H N {V N (t),u,y) = n + (V N (t),u)5 yu+ i_ + /j,-(V N (t),u)S ytU _i r 

where ^ + {V,u) = M/3( ^" ( M ^ (y) and ^(V,u) = ^y^Zl^v) ■ This 
jump law is also time-dependent through V/v(t). 

The deterministic system (D)takes the form: 
v(t) = f(v(t),g(t)) 

g{t) = {l-g{t))a{v{t))-g{t)p{v{t)) 
v(0) = v Q ; .g(0) = u 

In the sequel, we will be interested in the asymptotic behavior of the stochastic 
hybrid models (Sn) under the limit fluid assumption. Let us now recall what 
this assumption means. Let (xn) be a sequence of homogeneous Markov jump 
processes with state spaces En C R k 1 intensities Ajv(x) and jump law /ijv(x, dy). 
Define the flow as Fn(x) = Ajy(x) J En (z — x)fiN(x,dz). The fluid limit occurs 
if the flow admits a limit and if the second order moment of the jump size 
converges to zero when N — > oo. Our stochastic hybrid model is in the fluid 
limit since the jumps are of size 1/N and the intensity is proportional to N. 
This stems from the fact that we are modeling proportions in a population of 
independent agents. However, this independence assumption is not necessary 
to satisfy the fluid limit. 

2.2 Law of large numbers for stochastic hybrid systems 

We give here the first result concerning the convergence of the stochastic hybrid 
model (Sn), which is a functional law of large numbers on finite time windows. 

Theorem 2.1 () Let e > 0, S > 0, T > 0. Let us assume that the functions 
otij and f are C 1 , and satisfy the following condition: 
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(HI): the solution v of (D) is bounded on [0,T] and for all N > 1 the solution 
Vj\r(t) of (Sjsr) is uniformly bounded in N on [0,T]. 

Let X mlt a given initial condition for (D) and X = (v, g) its solution. Then 
there exists an initial condition for (5/v) o-nd No > such that V-/V > No, 



the solution X 



N 



(V N ,e 



(i) 



) satisfies, for all 1 < j < q and 1 < k < rj : 



sup \\V N (t)-v(t)\\ >6 

0<t<T 



< e 



sup |{e#h(t)-g£»(t)|>5 

0<t<T 



< e 



Moreover, if we define 
P N (T,A) = P 



sup J\V N (t) - v(t)\\ 2 + \\4\t) ~ S {j) (t)\\ 2 > A 



0<t<T 



there exist two constants B(T) > and C > such that for A sufficiently small: 

1 Ae- S ( T ) T 
lim sup - log Pn(T, A) < — — (2) 

Moreover if 



(H2): assumption (HI) holds true on [0, +oo[ 
then the constant B(T) — BT is proportional to T . 

Interpretation of the convergence speed We have obtained in (2.1) an 
upper bound for the convergence speed which can help to answer the following 
issue. Given a number of channels TV, given an error A and a confidence proba- 
bility 1 — p (e.g p — 0.01), the time window [0, T] for which we can be sure (up to 
probability 1—p) that the distance between the stochastic and the deterministic 
solutions (starting at the same point) is less than A is given by (2.1). In section 
3.3, we show numerical simulation results illustrating the obtained bound for 
the convergence speed for the stochastic Hodgkin- Huxley model. 

Remark Assumption (H2) and thus {HI), are satisfied for most neuron mod- 
els, for instance for the Hodgkin- Huxley (HH) model 0]. 

2.3 Langevin approximation 

Our second result is a central limit theorem that provides a way to build a 
diffusion or Langevin approximation of the solution of the stochastic hybrid 
system (S N ). For X = (v, e) e R p x R d , for 1 < j < q, 1 < i, k < r h let 
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l<i<rj ,i^k 

flg=aS(X)e«+ag(X)e« 
X j!k (X)= £ Hg- 

l<z<rj 

As before, Xjv(t) = (VAr(i), ejv(i)) € R p X R d is the solution of the stochastic 
hybrid model (Sjv). 

Let ifo(i) = {(iZ$) k (t)} i<i<«, i<fc<rj with Rffi e R rj be defined as 
(i^ } ) fc (t) = ViV ({e#} fc (t) - {e$} fc (0) - jf b jtk (X N (s))ds 

Theorem 2.2 Under the same hypotheses as in Theorem 2.1, the process Rm 
converges in law, as N — > oo, to the process R = {(R^')k(t)}i<j<q, i<fe<rj with: 

R {]) (t) = [ a (j \X(s))dWi 

where 



• X = (v, g) is £/ie solution of the deterministic model (D) with initial 
condition X lnlt = X™ u = X , 

• W J are independent standard r.j- dimensional Brownian motions, 

• (X) is the square root of matrix G^{X) s.t., for 1 < k,l < ry. 

1 G k J }(X) = Hij(X)=G^ k (X), l?k 

This theorem leads to the following degenerate diffusion approximation Xn = 
(Vn,Sn) € R p x R d , for TV sufficiently large: 



, 99 , i dV N = f(X N (t))dt 

[Z - Z) i dg« =b j (X N (t))dt+ -j=a^(X N (t))dW? 



where bj(X) is the vector {b^k{X))x<k< rj £ R rj , 1 < J < g. 

Note that this approximation may not have the same properties as the original 

process, even in the limit N — > oo (when considering for instance large deviations 

US)- 
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2.4 Functional central limit theorem and exit problem 

Let Xn = (Viv,e^, be the solution of the stochastic model (Sn) and 

X = (v, gW, .., g^) the solution of the deterministic system (D) with identical 
initial condition X init = X™ u = Xq G RP+ d . 
Consider the (p + d)-dimensional processes: 



( & \ 



N 



( v N 
J 1 ) 



-N 



V \ 

p-(l) 



With this setting, we have the following result: 

Theorem 2.3 Under the same hypotheses as in Theorem 2.1 the process Zn 
converges in law, as N — »■ oo fo i/ie process 



Z = 



( Y \ 

p(i) 

\ p(«) / 



whose characteristic function *&(t,6) = E [e t<e ' z W > ] satisfies the following 
equation: 



U) db itk 9* _ 1 ^ (J) (J) (3) U V V(9 ^ 3 * 



% -EjEEtf 

j=i I ;sl fe=i 



m=l igL 



and 3 g^ fc are evaluated at X(t), for 1 < m < p, 1 < j < 



a, 1 < k < rj and I G L, wif/i = ((0 m )i< m < p , (#£ )i<7<g, i<k<r 4 ) = {&l)leL, 
and L = {(m)i< m < p , (j, fc)i<j<?, i<fc<r 3 }- 

Let us define Zjy := (Yat, P/v) := v/iV (l^v - G R p x R d , where we recall 

that Ajv is the Langevin approximation defined in (2.2). Thus, for 1 < j < q, 
and 1 < k < rj-. 

dY N = 



dP 



pj.k 



N 



N(f(X N ) - f(X))dt 
N{b 3 , k {X N ) - b jlk (X))dt + a^(X N )dWi 



As an asymptotic linearization of Zm, we define the diffusion process 
(7, 7r) G R p x R d by: 



d7„ 



^56 



fc'=l 
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Theorem 2.4 The processes Z and 6 have the same law. 

A computation of the moments equations for the limit process Z for the example 
of section 2.1 is provided in Appendix B. 

Using the central limit theorem 2.3, we provide in the next theorem a characteri- 
zation of the fluctuations of first-exit time and location for the stochastic hybrid 
process X?q. Let <f> : W +d — > W +d be continuously differentiable. Define 

t n := inf{i > 0; <f>(X N (t)) < 0} 

t := inf{i > 0;4>(X(t)) < 0} 

the first passage times through <f> = respectively for the stochastic hybrid 
process and for its deterministic limit. 

Theorem 2.5 Assume the initial condition X(0) satisfies <j>(X(0)) > 0. Sup- 
pose t < oo and V cf>(X (t)) .F (X (t)) < 0. Denote the random variable 

V^X(t)).Z(t) 
1 > • V<A(X(r)).F(X(r)) 

Then the following convergences in law hold when N — > oo : 

VN(tn - t) -> 7r(r) 
^(^(tat) - X(r)) -> Z(t) + 7r(r)F(A(r)) 

3 Application to neuron models 

In this section, we show how our previous theorems can be applied to standard 
neuron models taking into account ion channel stochasticity. 

3.1 Kinetic formalism in neuron modelling 

Kinetic models can be used in many parts of nervous system modelling, such 
as in ion channel kinetics, synapse and neurotransmitters release modelling. 
Deterministic kinetic equations are obtained as a limit of discrete stochastic 
models (hybrid or not) as the population size, often the number of channels, is 
large. 

As already mentioned in the introduction, compared to our general formal- 
ism, the stochastic individuals are the ion channels and the global variable Vn 
the voltage potential. Constituted of several subunits called gates, voltage-gated 
ion channels are metastable molecular devices that can open and close. There 
exist different types of channels according to the kind of ions, and they are 
distributed within the neuron membrane (soma, axon hillock, nodes of Ranvier, 
dendritic spines) with heterogeneous densities. 

In what follows, we consider the model of Hodgkin and Huxley, which has 
been extended in different ways to include stochastic ion channels. In numerical 
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studies, different versions have been used, from a two-state gating interpretation 
e.g. [31] to a multistate Markov scheme [7J [28]. In [32], two of these models 
are compared, one with a complete multistate Markov scheme, and the other 
inspired from [53] with a multistate scheme for the sodium ion and a two-state 
gating for the potassium ion. We are here considering only single-compartment 
neuron, but in order to deal with spatial heterogeneities of axonal or ion channels 
properties for instance, multi-compartmental models can be introduced as a 
discretized description of the spatial neuron, with Ohm's Law coupling between 
compartments. 

3.2 Application of the law of large numbers to Hodgkin 
and Huxley model 

Classically, the Hodgkin-Huxley model (HH) is the set of non-linear differential 
equations (3.1-3.4) which was introduced in [TB] to explain the ionic mechanisms 
behind action potentials in the squid giant axon. 

/ - 9l(V - V L ) - g Na m 3 h(V - V Na ) ~ g K n A {V - V K ) (3) 

{1 - m)a m (V) - mp m {V) (4) 

{1 - h)a h {V) - hp h {V) (5) 

(1 - n)a n (V) - np n {V) (6) 

where / is the input current, C m — IfiF/cm 2 is the a capacitance corresponding 
to the lipid bilayer of the membrane, ql — 0.3mS/cm 2 7 gNa = 120mS / cm 2 , 
gK = 36ms/ cm 2 are maximum conductances and Vl = 10.6mV, VNa — llbmV, 
Vk — — 12mV are resting potentials, respectively for the leak, sodium and 
potassium currents. The functions a x , f3 x for x = m,n,h are opening and 
closing rates for the voltage-gated ion channels (see [TB] for their expression). 
The dynamics of this dynamical system can be very complex as shown in |14j . 
but for our purpose let us describe only schematically the behavior of this system 
as the parameter / is varied (see [37] for more details). First for all I there exists 
a unique equilibrium point. For < I < I\ 9.8/iA/cm 2 , this equilibrium point 
is stable, and for Iq < I < l\ where Jo « 6.3fj,A/cm 2 this equilibrium coexists 
with a stable limit cycle and possibly many unstable limit cycles. At I = I\ 
and I — h occur two Hopf bifurcations. For I\ < I < I2 « 153[iA/cm 2 , 
the equilibrium point is unstable and coexists with a stable limit cycle. For 
I > I2, there are no more limit cycles, and the equilibrium point is stable. 
This bifurcation structure can be roughly interpreted as follows : for / small 
the system converges to an equilibrium point, and for I sufficiently large, the 
system admits a large amplitude periodic solution, corresponding to an infinite 
sequence of action potentials or spikes and the spiking frequency is modulated 
by the input current /. 



C 



dV 
l ~dt 
dm 

~~dt 
<lh_ 

~dt 

dn 

~dt 



12 



There are two stochastic interpretations of the Hodgkin- Huxley model involving 
either a multistate Markov model or a two state gating model. We now present 
them briefly and we apply our theorems to each of them. A comparison of 
these deterministic limits obtained for these models is given in Appendix A and 
establishes an equivalence between the deterministic versions as soon as initial 
conditions satisfy a combinatorial coincidence relationship. This question is 
further studied in [19], where the reduction of the law of jump Markov processes 
to invariant manifolds is investigated. 



Multistate Markov model This model has two types of ion channels : one 
for sodium and the other for potassium. The kinetic scheme describing the 
Markov jump process for one single potassium channel is the following: 

©^©^©^©^© 



k ^ ' 2£>, 
And for sodium channel 



4 P, 




All the coefficients in these two schemes are actually functions of the membrane 
potential, and can be found in [16]. The state spaces are: 

Ex = {no,ni,n2,n 3 ,n 4 } 

E 2 = {mohi,mihi,m2hi,m3hi,moho,miho,m2ho,m3ho} 

The open states are respectively n.4 (r\ = 5) and m^hi (r 2 = 8). The proportion 
of open potassium channels is denoted by :— {e^'}„ 4 and the proportion 

(2.) (2.) 

of open sodium channels by u y N ' := {e K N } m3 h ± - In this model, The membrane 
potential dynamic is given by the equation: 

V N (t) = -g N auf(t)(V N {t) - V Na )-g K u%\t){V N {t) - V K )-g L {V(t) -V L ) + I 

where I € R is a constant applied current. With the notations of the previous 
sections, f(v, u«, u^) = -g Na u^ ' (v-V Na )-g K uW ' {V N (t)-V K )-g L (v-V L ) + 
I and for example, cx^Av) = 3a n (v) if k = n%, j — n 2 and oi^\(v) — 2j3 rn {v) if 



k = m 2 ho, j — miho. 

Applying Theorem 2.1, we obtain a deterministic version of the stochastic 
Hodgkin- Huxley model when N — > 00, provided we choose the initial conditions 
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appropriately: 

' v = ~9Nae^(t)(v(t) - V Na )-g K eW(t)(v(t) - V K ) - g L (V(t) - V L ) + 1 

^(0) = V init 

{ g «(o) = gE ) li 

with — gif , and where the rate functions a~m. p are given in the above 
schemes, for j £ {1,2}. 

Two-state gating model Another way of building a stochastic Hodgkin- 
Huxley model is to consider that the channels can be decomposed into indepen- 
dent gates. Each gate can be either open (state 1) or closed (state 0): 

I) 1 

with z € {to, n, h}. The channel is open when all gates are open. 
So, here, q = 3 and E\ = E 2 = E a = {0, 1}. If we denote u^\t) :— {e^}i the 
proportion of open gates z, for z G {to, n, h}, the membrane potential dynamic 
is then given by: 

V N (t) = -9Na(u ( ™\t)) 3 u%\t)(V N (t)-V Na ) 

- g K (u%\t))\V N (t) - V Na ) - g L (V N (t) - V L ) + I 

which corresponds to /(«,«("•) = ~gNa{u (m) fu^\v-V Na )-g K {u {n) ) A {v~ 

Vno) — 3l (v — Vl) + 1- In Figure 1, we give a sample trajectory of this two-state 
gating stochastic Hodgkin-Huxley system. 

Applying Theorem 2.1 gives the classical formulation of the 4-dimensional 
Hodgkin-Huxley model: 

v = - 9Na (uW(t)) 3 uW(t)(v(t) - V Na ) - g K (u^{t))\v(t) - V Na ) 
-g L (V{t)-V L )+I 

ii {z \t) = (I - (t))a z (v(t)) ~ (t)/3 z {v(t)), z e {m,n,h} 
3.3 Exponential convergence speed 

We illustrate by numerical simulations the upper bound obtained in (2.1) for 
the stochastic Hodgkin-Huxley model with a two-state gating scheme. The 
number of sodium channels N^ a and potassium channels Nk are proportional 
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5 10 15 20 25 30 35 40 45 




5 10 15 20 25 30 35 40 45 50 




Figure 1: Illustrative sample trajectory of a two-state gating stochastic Hodgkin- 
Huxley model with N = 20 (cf. section 3.). Top: variables m, h for the sodium 
channel (without unit). Middle: variable n for potassium channel (without 
unit). Down: variable V for membrane potential (unit: mV). Abscissa: time 
(arbitrary units). Since m, n and h correspond to proportions of open gates, if 
one of them is equal to 1, it means that all the corresponding gates are open. 
An increase in the membrane potential V causes an increase in the proportion 
of open m (sodium) gates, which in turn implies an increase of V. This positive 
feedback results in a spike initiation. Meanwhile, a further increase of V is 
followed by a decrease of the deactivation variable h, which closes the sodium 
channels. This inhibition effect acts at a slower time-scale, enabling a decrease 
of V. This decrease is strengthened by the dynamic of variable n (proportion 
of open potassium gates) . 

to the area S of the membrane patch. Thus, instead of N, S will denote the 
size parameter. For the squid giant axon, the estimated densities for the ion 
channels used in the simulations are pNa = 60^tm~ 2 and pm k — 18fim^ 2 . 
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We now display the results of numerical simulations of 




using Monte-Carlo simulations. We recall that from (2.1): 

1 Ae- BT2 
limsup - log P S (T, A) < — - = C{T) 

S->oo O 

In Fig. 2, the simulation estimations of Cs{T) = i logPg(T, A) are shown for 
different values of T and S and can be compared to the theoretical bound C(T). 
Simulations are made without input current, meaning that the stochastic solu- 
tion is supposed to fluctuate around the equilibrium point of the deterministic 
system in a neighborhood of size proportional to S^ 1 ^ 2 . When S increases, the 
simulation curve Cs(T) is expected to pass below the theoretical bound C(T). 
For higher input currents, still subthreshold (I < I e ), but close to the bifurca- 
tion, channel noise will induce spontaneous action potentials. For appropriate 
A, the probability Ps(T, A) can be interpreted as the probability that the first 
spontaneous action potential (SAP) occurs before time T. Thus the convergence 
speed bound gives an upper bound of the repartition function of this first SAP 
time. 

For higher input currents I > I c , the deterministic solution will be attracted 
by a stable limit cycle, which corresponds to repetitive action potentials. In 
this case, channel noise can introduce a jitter in the spiking times. Thus, if one 
considers the supremum of the errors between the stochastic and the determin- 
istic solutions, this supremum will be quite large (approximately the size of an 
action potential) as soon as the difference between the stochastic spiking times 
and the deterministic ones is of order the time course of an action potential (2 
ms). Thus, the supremum of the difference is not appropriate here and we will 
see in the following section how to quantify the impact of channel noise on the 
spiking frequency. 

3.4 Application of the central limit theorems 

In this section, we show how to investigate the fluctuations around a stable fixed 
point (sub-threshold fluctuations) and the fluctuations around a stable limit 
cycle (firing rate fluctuations) using Theorem 2.3. Let us consider a class of two- 
dimensional models, corresponding to the Example of section 2.1. This class 
contains reductions of the previous two-state gating Hodgkin- Huxley model, or 
other models such as the Morris-Lecar model [25]. Consider the process, with 
the notations of the Example: 

( Y N \ = ( VN(V N -V) \ 
\Pn J \ VN(e N -g) ) 
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-4.0e-03H 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 

50 100 150 200 250 300 350 400 45C 

T 

Figure 2: Simulation results of the Hodgkin-Huxley model with a "two-state 
gating" scheme with input current 1 = 0: this figure shows the quantity 
g In Ps(T, A) as a function of T, where S is the area of the patch, and is thus 
proportional to N. Stars : S = 250^m 2 (corresponding to N Na = 15000 and 
N K = 5000) ; Empty boxes : S = 500/xm 2 ; Crosses: S = 750/um 2 . Lines are 
guide for the eye. 



with initial conditions (Pjv(0), Yjv(0)) = (0,0). Then the 2-dimensional process 
Zn = (Pn,Yn) converges in law, as TV — > oo, towards the process Z = (P,Y), 
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whose characteristic function is given by: 



E 



e i(etP(t)+e 2 Y(t)) 



A' C /2 s * 

Thus, defining E s the square root matrix of T s := ( nl s , SJ, I , for < s < t, 



Z can be written as a gaussian diffusion process: 



Z t = / Z s dW s 



CJ2 B> 



o 

where W is a standard two-dimensional brownian motion 0. 
From the equation for the characteristic function obtained in Theorem 2.3, 
one derives that the triple y = (A, B,C) is solution of the system y t = M t yt + E t 
defined as: 

-\\{V,u) 

| (M) 








(V 


H 


V c t ) 







o K " 







2/; fu 


Is 


V'u 







with initial conditions (0,0,0), and X(v,u) = — u)a(v) + u/3(v). The par- 
tial derivatives f' v , f' u , b' v , b' u and A are evaluated at the deterministic solution 
(Vufjt) ■ 

We remark that, if J is the Jacobian matrix at the point (Vt,gt), and if its 
spectrum is sp(J) = {Ai, A 2 } then the spectrum of M is Sp(M) = {2Ai, 2A 2 , Ai + 
A2}. Two different situations can be considered: 

• Starting from a fixed point (Vo, uq) of the deterministic system, the matrix 
M t = M(Vt,u t ) and the vector E t = E(V t ,u t ) are constant. One can 
derive an explicit analytical solution diagonalizing the matrix M. The 
time evolution for the variance and covariance of the difference between 
the deterministic solution and the stochastic one then depends on the 
stability (Ai, A2) of the considered fixed point. 

• Around a stable limit cycle (periodic firing): M t and E t are T-periodic 
functions. Using suitable coordinates and following Floquet's theory (see 
[3]), stability would be given by the spectrum of the solver R(T) : (A , B , Co) 
(At, Bt, Ct)- As explained in [17], even if the real parts of the eigenval- 
ues of the jacobian matrix are strictly negative for all time, unstable solu- 
tions may exist. In section 3.5 we investigate numerically the fluctuations 
around a stable limit cycle for the Morris-Lecar system. 



4 The condition that the matrix T s admits a real square root matrix can be reduced to 
K + B' s < because one can show that det(r s ) = A' S B' S - C' a 2 /4 = for all s > 0. This 
condition is thus always satisfied because : A' + B' < 0, A' s and B' s have the same sign, and 
(A' a , B' a , C" s ) cannot cross (0,0,0) by uniqueness of the solution of z' = Mz (satisfied by y'). 
The computation of the matrix E s gives: 

_ ^-2(A' S + B' S ) 
A> + B> 
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If we consider 

/ Y N \ = ( VN (V N - V) \ 
V Pn J \ </N(u N -u) j 

where (Vat, mat) is the Langevin approximation, then the moments equations, 
written for the linearized version around the deterministic solution, give the 
same matrix T s at the limit N — > oo. But for finite N the linearized process is 
not gaussian (see Appendix B). Thus, our mathematical result can be directly 
related to the simulations results obtained in [32:: in this paper simulations of 
two neuron models with a large number of stochastic ion channels are made, 
and the fluctuations of the membrane potential below threshold exhibit approx- 
imately gaussian distributions, but only for a certain range of resting potentials. 
For smaller resting potentials, the shape of the distribution remained unclear 
as it was more difficult to compute. Our approach shows that, at finite N, for 
any range of the resting potentials the distribution is non-gaussian, but when 
N — > oo, the distribution tends to a gaussian, which corresponds to the approx- 
imate gaussian distribution observed in the simulations of [32J. 

3.5 Quantifying the effect of channel noise on neural cod- 
ing 

Neurons encode incoming signals into trains of stereotyped pulses referred to as 
action potentials (APs). It is the mean firing frequency, that is the number of 
APs within a given time window, and the timing of the APs that are the main 
conveyors of information in nervous systems. Channel noise due to the seemingly 
random fluctuations in the opening and closing times of transmembranar ion 
channels induces jitter in the AP timing and consequently in the mean firing 
frequency as well. We show in the next subsections how our results can be 
applied to quantify these phenomena. The impact of channel noise on frequency 
coding is investigated in sec 3.5.1 and on spike timing coding in section 3.5.2. We 
close this section by some remarks concerning non-markovian processes arising 
when considering synaptic transmission in sec. 3. 5. 3. 

3.5.1 Numerical study of the variance of spiking rate for Morris- 
Lecar model 

In this subsection, applying Theorem 2.3 to the Morris-Lecar system, we in- 
vestigate the impact of channel noise on the variance of the firing frequency. 
The Morris-Lecar system was introduced in [25] to account for various oscillat- 
ing states in the barnacle giant muscle fiber. We denote by X = (V, m, n) the 
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solution of: 



cj^- = I-gL(V-V L )-gcam(V-V C a)~gKn(V-V K ):=F v (Xl7) 

^ = A ro (^)(M 00 (^) - m) := F m (X) (8) 

— = KiV^NUV) - n) :=F n (X) (9) 

whrer X m {V) = cosh((V-V 1 )/2V 2 ), A n (V r ) = </>„*cosh((V- V 3 )/2^), Af oc (V r ) = 
(1 + tanh[(V - 7i)/V 2 )])/2 and ^(V) = (1 + tanh[(V - V 3 )/V 4 )})/2. We 
introduce as in the previous sections a stochastic version Xn of this model 
with stochastic ion channels, replacing the differential equation for m and n by 
birth-and-death processes with voltage-dependent opening rates a m — A m Moo, 
ot n = AniVoo and closing rates j3 n = A n (l — A^). According to the parameters of 
the model, the deterministic system (3.5 — 3.7) may have a stable limit cycle x LC 
for some values of I €E [I m in, Imax] (see |25j). This corresponds to a phenomenon 
of regular spiking, characterized by its rate. Assuming that the time length of 
a spike is almost constant, we suggest a proxy for this spiking rate: 

1 f T 

r ( T ) ■= 7p J o <j>th{x{s))ds 

where <f>th is a sigmoid threshold function. In a similar way, we define the 
stochastic spiking rate by: 



r 



r N (T) := i / <f> th (X N (s))ds 



As a candidate for 4> t h, we choose 4>th(V) := - a(v Jv\i where c and Vth are 
two parameters. 

A consequence of the central limit theorem for Xn is the following weak con- 
vergence: 

1 ' T 



TV [r N (T) - r(T)] i?(T) - - / Z(*).V& h (a;(«))ds 



where Z is the weak limit of yN [Xn — x]: 

Z(s) = ( Z{u)dW u 



R(T) is a gaussian random variable with zero mean. For simplicity we consider 
the case where <pth is only a function of the membrane potential V . Then the 
variance of R(T) is: 



4(T) = E [R(Tf] =^J^ J's v ^W^V^))dJ^{y{s))da (10) 
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where S v (s) — £1,1 (s) is the variance of y/N(Vff(s) — V(s)). 
To estimate numerically the variance er|j(T), the first step is to determine nu- 
merically the limit cycle, then solve the moment equations (Appendix C) and 
immediately deduce E(s). Thus the variance a 2 R can be computed using formula 
(3.8) without any stochastic simulation. In Fig. [3] we show our numerical re- 
sults, where we plot in C-F., as a function of the input current /, the normalized 
variance £(T) defined as: 




D Input current I E Input current I F Input current I 

Figure 3: Impact of channel noise on the spiking rate. First row (ABC) : Class 
I regime. Second row (DEF) Class II regime. [A-D]. Deterministic rate r(T) 
versus input current I (p,A/ cm 2 ). [B-E]. Variance <j\(T) versus input current / 
(fiA/cm 2 ). [C-F]. Normalized variance £(T) versus input current I (/j,A/cm 2 ). 
Parameters : for all figures T = 2000ms, c = 10, V th = OmV, C m = 20fiF/cm 2 , 
Vi = OmV, V 2 = 15mV, V 3 = WmV, gca = AmS/cm 2 , g K = 8mS/cm 2 , 
g L = 2mS/cm 2 , V K = -70mV, V L = -50mF, V Ca = lOOmV, tf> n = 0.1 and 
V 4 = 20mV for Class II (DEF) and V 4 = WmV for Class I (ABC). 

Comments The value of £ (T) depends on a combination of the linear stability 
along the cycle and on the variance of the noise (which is multiplicative) along 
the cycle. If one wants to have the quantity E [( rjv (^~J'( T )) 1 f order 1, then 
the number N of channels should be of order £(T). Interestingly, this gives 
much smaller values for Class II than for Class I regime. In both cases, it 
corresponds to a reasonably small number of channels when / is not too close 
from bifurcation points. 
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3.5.2 Impact of channel noise on latency coding in the Morris-Lecar 
model 



Whereas frequency coding requires an integration of the input signal over a 
relatively long time, individual spike time coding does not require such an in- 
tegration. The time to first spike, called latency, depends on the value of the 
suprathrcshold input. Thus it may have an interpretation in term of neural cod- 
ing, and it has been shown in several sensory systems |36j that the first spike 
latency carries information. For example, a recent study [13] concerning the 
visual system suggests that it allows the retina to transfer rapidly new spatial 
information. Impact of external noise on latency coding have been investigated 
in numerical studies [9] with stochastic simulations. We apply Theorem 2.5 to 
the Morris-Lecar model to investigate the impact of internal channel noise on 
first spike time. We chose the parameters (see [3]) to obtain a Class I neuron 
model in the excitable regime. In this setting, there exists a unique steady 
state X* = (V* ,m* ,n*). Starting from this equilibrium point, the impact of 
an input at t = to is equivalent to an instantaneous shift of the membrane po- 
tential V* -> V* + A, where A > is the amplitude of this shift. Eventually 
the system goes back to its steady state, but if A is higher than a threshold 
A t h then a spike is emitted before going back to the steady state, whereas if 
A is lower than A t h no spike is emitted. For A > A t h, we define the latency 
time T(A) as the elapsed time between to and the spike. More precisely, let 
(Va(£), 771,4 (t), for t > t be the solution of Morris-Lecar equations with 

initial conditions X(t ) = (V(t ) — V* + A,m(t Q ) = m*,n(to) = n*). We de- 
fine a spike as a passage of the membrane potential Vk(i) through a threshold 
Vth- Then, with to = for simplicity, the latency time T{A) can be written 
as T(A) := m£{t > 0; V A (t) > V th }- As shown in FiggjA, the more A > A th 
is close to A t h, the longer is the latency time T(A). The same setting can be 
extended in the stochastic case, defining a random variable Tm{A). Applying 
Theorem 2.5, with 4>(V, m, n) — Vth — V, we express the variance P(A) of the 
limit of y/N(T N (A) - T(A)) as N -> oo: 

_ S V {T(A)) 
P{A) ~ F V {X(T{AW [ ] 

In (3.9), S V (T(A)) is the variance of the ^-component Z v of Z, where we recall 
that Z is the limit of y/N(X N -X) (see Theorem 2.3). The value of S V (T(A)) is 
obtained from the numerical integration of the moments equations (..). The re- 
sults are displayed in Figf?l where the variance P(A) and a normalized variance 
P(A)/T(A) 2 are plotted against the amplitude A (g|B). In SJD the variance 
P{A) is plotted against the latency time T(A) (jUD). From (3.9), it appears 
that P(A) is determined by two distinct contributions : the variance S V (T(A)) 
(g]E) and the crossing speed F(X(T(A))) (j3]F) which actually does not influ- 
ence much the variance P(A). One way to interpret the results is the following: 
if N is large, of order P(A), then E[(T N (A) - T(A)) 2 ] is of order 1. Thus, as 
an illustration, in order to keep E[(Tn(A) — T(A)) 2 ] of order 1, the required 
number of channels would be of order 10 2 for a latency time of 10ms and of 
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Amplitude A (mV) Amplitude A (mV) 

Figure 4: Impact of channel noise on latency coding. [A]. Latency time T(A) 
versus amplitude A. [B]. Variance P{A) versus amplitude A. [C]. Normal- 
ized variance P(A)/T(A) 2 versus amplitude A. [D]. Variance P(A) versus la- 
tency time T(A). [E]. Variance S V (A) versus amplitude A. [F]. Crossing speed 
F V (X(T(A))) versus amplitude A. Same parameters as in [3] Class I, with input 
current I — 32/xA/cm 2 . 

order 10 5 for latency time of 60ms. 

3.5.3 Synaptic transmission and non-markovian processes 

In section 3.5.1, the quantity of interest was the firing frequency. However, 
the synaptic transmission between a neuron A and a neuron B has its own 
time scales. Therefore, neuron I?'s input, called post-synaptic potential 1 i> A ~ >B , 
may be modeled as a functional of neuron A's membrane potential {VA(t)}t>a- 
Although synaptic transmission is presumably a non-linear process, one can 
consider as a first approximation (cf. [21]) that the process of interest is obtained 
directly by the convolution of the process Va with some kernel K A ^ B : 

^ B {t)= [ K A ^ B (t- s)V(s)ds 
Jo 

The mathematical analysis of the impact of channel noise on this variable 
can be done in the light of theorems 2.1 and 2.3. Using the general nota- 
tions for the stochastic process and its deterministic limit, we define \I/./v(t) = 
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J*K(t,s)X N (s)ds and #(*) = /„* K(t, s)x(s)ds. 

Law of large numbers Define 

S N (T) = sup \y N (t)-V(t)\ 2 
se[o,T] 

Clearly, using Cauchy-Schwartz inequality: 

P(S N (T) > A) < P^T^T^A) 

with 

rj(T) = T sup / s)| 2 ds 
te[o,T] Jo 

The convergence of ^jv to \t with the same kind of exponential convergence 
speed is thus a direct consequence of Theorem 2.1. 

Gaussian fluctuations We know (Theorem 2.3) that \ /r N(Xff — x) converges 
weakly to the diffusion H(t) = J Q R(u)dW u . As a consequence, fi^ = V^V (^jv — 
'J/) converges also weakly, to the following process: 

fi(t) = ^ if(t,s) ^ i?(u)dW M ^ 
With an integration by part, one can rewrite: 

n(t) = / z(t, s)dw s 

Jo 

with 

Z(t, s) = / K(t,u)duR(s) 

J s 

The process fi is gaussian and one can easily compute its variance as J Q Z(t, s) 2 ds. 
However, it is non markovian, and some issues concerning the first hitting times 
of such processes are solved in [33J . 



4 Proof of the law of large numbers 

In this section we give the proof for Theorem 2.1. This proof is inspired 
from pQ, except for the exponential martingale bound. In order to simplify the 
notation and to make the arguments clearer and more intuitive, we write the 
proof for the case of a single channel type with state space {0, 1} and transition 
rates given by the scheme: 

In this case, the stochastic model (S%) is: 

V N (t) = f(V N (t),u N {t));V N (0) = V 
1 N 

fc=l 
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«(v N > 



where ?4 : — > 1 with rate a(Viv(t)) and 1 — > with rate /3(V/v(i)), for all 
1 < k < N 

The deterministic solution (v, u) satisfies: 



In order to complete the proof, few slight changes in the notation can be done: 

• in order to work with more general jump Markov processes with finite state 
space, essentially all the expressions of the form So(u)a(v) — Si(u)a(v) 
should be replaced by 



• in order to include q different channel types (different ions), one should 
just write the same arguments for all the q processes {e^\t)} for 1 < j < q 

and include all the ||e^(i) — e^(f)|| for 1 < j < q in the function f(t) of 
Gronwall lemma application in section 3.4. 

4.1 Decomposition in a martingale part and a finite vari- 
ation part 

Decomposition We decompose the difference between the stochastic and the 
deterministic processes as a sum of a martingale part Mm and a finite variation 
part Qn as follows: 



v(t) = /(«(*), u(t)) 

u(t) = (1 - u(t))a(v(t)) - u(t)P{v(t)) 
v(0) = v ; u(0) = u 




[u N (t) - ujv(O)] - [u(t) - u(0)] = M N (t) + / Q N (s)ds 




where we define: 



Ow(t) = E [Wu^MVkKt)) - M^MM*))] - u{t) 



M N (t) = [u N {t) - u{t)] - [u N {0) - u(0)] 
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Lemma As defined above, (Mjv(t)) is a {i^j-martingale. 

Proof For h > 0, define AM N (t, h) = £E [M w (t + h) - M N (t)\F t ], then: 



AM N (t,h) 



1 1 

ft iV 



N 

E E N u 



JV 



-E 



*i(t4°)i^ 



/■t+h I ^ 

L »=i 



ft 



[u(t + ft) — u 



1 /■*+/> 



u(s)ds 



The last line converges clearly to as ft — > 0, and the two first terms compensate 
as ft — > 0. So we have: 



Therefore 



lim ^E [Mjv(t + ft) - Mjv(t)|F t ] = 

h-yO a 



-E[M N (t + s)\F t ]\ s=o =0 



By dominated convergence we have: 
^-E[M N (t + s)\F t ]\ s=S0 =E 



— ~E[M t+S0+u \F t+SQ ]\ u = \F t 



Finally: 



E[M N (t + h)\F t ] = Cst = M N (t) 



4.2 Martingale bound 

In this part we want to obtain a bound in probability for the martingale part. 
We introduce the jump measure and the associated compensator: 
We define two random measures on ]0, T] x {0,1}: 



E 



• jump measure : Kt = > 5,. (*), 

te]o,T],u^V«f2 

• compensator : 

Vi{dt,dy) = b{V N {t))5 1 {uf)5 {y) + a{V N {t))5 {uf)5 1 {y) 



dt 



We can rewrite Qn(s) and Mjy(t): 

I Q N (s)ds = V" / (^i(y) - ^(^'.^^(ds,^) - / u(s)ds 

JO iV ~^J]l),T]x{0,l} JO 

= / (Ml/) - M«t-))(Ki - 

~lJ]0,T]x{0,l} 



Then we have the following proposition: 
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Proposition Let T > 0, e > 0, 8 > 0. Then there exists Nq such that 



sup M N (t) 2 > 8 

0<t<T 



< e 



Proof Let us first recall that from standard results about residual processes 
([18]) we have: 

i N r r 

E[M N (t) 2 ] = / {8 1 {y)-8 l {uf)f{K i -y i ){ds,dy) 

iV j=i L^]0,T]x{0,l} 



1 ^ 



0,T] 



/3(^v( S ))5 1 (u«)+a(^ v ( S )),5o(u«)d S 



Therefore, we can get a bound for E [Mjv(i) 2 ] : 

E[M N (t) 2 ] ^d^maxdlalU.lljSIU) 

where ||a||oo and ||/3||oo are finite because a and /3 are continuous and assump- 
tion (HI). We then use Chebychev inequality and Doob inequality for L 2 mar- 
tingales: 



sup M N (t) 2 > 8 

0<t<T 



4- 



sup M N (t) 2 

0<t<T 



< -E [M N {tf 



and E [M N (t) 2 ] < f for all N > N .f< 



In order to obtain a better estimate for the convergence rate, we derive 
here an exponential bound for the martingale part. Our proof is inspired from 
techniques developed in [5]. 



Proposition Let T > 0,?/ > 0. There exists a constant C v such that for all 
8 e]0, r)C n T[: 

8 2 N 



sup \M N (t)\ > 8 

0<t<T 



< 2 exp 



2C n T 



Proof We define, for x — (u,v), 9 e R: 

m N (x,9) = [ e 6v \ N (x)n N (x,dy) = N\(x){e 9/N n + (x) + e- e / N n-(x)} 



4>n{x,9) = I ' [e 0v -l-9y]X N {x) m {x,dy) 

JR 
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The second equality stems from integration by part. And, if \6\ < Nrj, 



d 2 mjy 



(x,r6)\ = 
1 Cr, a 2 



1 96» 2 

So, \(p N {x,6)\ < hlT 02 - Let us define 



< 



~W 



Z e N (t) = cxp[eMjv(t) - f <j> N ((u N {s),V N {s)),e)ds] 
Jo 

{Z e N {t)) is a martingale thanks to Doleans Formula: 

Z%(t) = l+ I f Z* N (s-)[e*y-l}(n-v)(ds,dy) 
Jo Jr 



Then we note r = mi{t;M N {t) > 6}. On {r < t},Z e N (r) > exp{c5e - 
And by optional stopping theorem: 

E[^(min(t,r))] = E[Z^(0)] > E[^(r)l T < t ] > P(r < t)exp{5e- ^} 



So, P 



sup I Mjv (t) I > <5 

0<*<T 



= P(r<T)<exp{-(5e+^}. 



Finally when 5 €]0, ryC^Tf, with e = and applying the same argument to 
— Mjv(t) we get the result.® 

4.3 Finite Variation Part 

In this section we use the Lispchitz property of a and (5 to provide a bound for 
the finite variation part, in order to apply later Gronwall Lemma. 

Lemma There exists C\ > independent of N such that: 

|Qjv(*))| < C (|ujv(i) - u{t)\ + \V N (t) - v{t)\) 

Proof 

1 N 

i=l 
1 ^ 
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Let us start with the second term of the difference, called Q 1 ^ : 
1 N 

1 ^ 

= ^ E M^WM*)) - + «(t) (/?(^v (t)) - ^(«(*))) 

i=l 

= f3(V N (t)) (u N (t) - u(t)) + u(t) (f3(V N (t)) - f3(v(t))) 

e[o,i] 

Then, 

IQ 1 ^ ! < \\f3\\°o\u N (t)-u(t)\+K \V-v(t)\ 

where Kp is the Lipschitz coefficient of (i. We do the same for the other term 
of the difference: 

IQ ^ 1 ! < IMUM^-uWI + tfalM*)-^)! 

So the proof is complete, with C\ = max(||a||oo, ||/3||co, K a , Kp) <g> 

If more general transition rates a(v, u) and 0(v, u) depend on v and u, 

one would need to replace ||a||oo and ||/3||oo, respectively by Halloo + and 
||/3| loo +Kf\ where K^\k^ are the Lipschitz coefficients associated with 
the second variable u. 

4.4 Proof of theorem 2.1 

Law of large numbers We want to apply Gronwall Lemma to the function: 

f(t) = \V N (t) - v(t)\ 2 + \u N (t) - u(t)\ 2 

From the previous section we have a good control on the martingale term and 
the following estimate: 

Corollary There exists C2 > independent of N such that: 

\u N {t) - u{t)\ 2 < 4[\u N (0)-u{0)\ 2 + C 2 T f \u N (s) - u(s)\ 2 ds 

Jo 

+ C 2 T f \V N (s)-v(s)\ 2 +M N (tf] 
Jo 

Proof As u N (t) - u(t) = ujv(O) - u(0) + M N {t) + /„' Q N {s)ds and (x + y + 
z) 2 < A{x 2 + y 2 + z 2 ), the result is a direct application of the previous lemma 
and of Cauchy-Schwarz inequality.® We need now to work on |Vjv(i) — v(t)\ 2 , 
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using hypothesis (HI) , with K\ = sup sup 

N se[0,T] 

df 

sup sup — (V N (s),u N (s)) 

N s£[0,T] OU 

Between the jumps, we have: 



df 

— {V N (s),u N {s)) 



and Kn = 



d_ 

dt 

Thus 



(\V N (t)-v(t)\ 2 ) = 2(f(V N (t),u N (t))-f(v(t),u(t)))(V N (t)-v(t)) 



\V N (t)-v(t)\ 2 = 2/ [f(V N (s),u N (s))-f(v(s),u(s))}(V N (s)-v(s))ds 
Jo 

+ \V N (Q)-v \ 2 

< \V N (0) -v \ 2 + 2K X f \V N {s) - v(s)\ 2 ds 
Jo 

+ 2K 2 I \u N (s) -u(s)\\V N (s) -v(s)\ds 
Jo 

\V N (0)-v \ 2 + 2K! f \V N ( S )-v( S )\ 2 ds 
Jo 

f \u N (s)-u{ S )\ 2 ds + K 2 f \V N (s)-v(s)\ 2 ds 
Jo Jo 



< 



+ K; 



where we used successively Cauchy- Schwartz inequality and ab < \{a 2 + b 2 ). 
Putting together this inequality with the Corollary we obtain: 



/(*) < A + B 



I f{s)da 
Jo 



where B = B(T) = max(2 J ft'i(T) + K2(T), C 2 T) does not depend on N and is 
linear w.r.t T if (H2) holds, and 

A=\u N {i))-uv\ 2 + \V N {i))-v \ 2 +K A sup M 2 

0<s<T 

If we control the initial conditions, then, with the control we have on the mar- 
tingale part, A can be chosen arbitrarily small (with high probability) and we 
can conclude with Gronwall Lemma. 



Exponential convergence speed If the initial conditions are the same for 
the stochastic and deterministic model, we actually have a exponentially fast 
convergence, thanks to the exponential bound for the martingale part: there 
exists a constant C m > such that: 



lim sup -|- log P 



sup \V N (t) - v(t)\ 2 + \u N (t) - u(t)\ 2 > A 

0<t<T 



< -- 



Ae -B(T)T 

2KAC m T 
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5 Proof of the central limit theorems 



As before, we write the proofs for the case of a single channel type with state 
space {0, 1} and transition rates given by the scheme: 




5.1 Langevin approximation 

In this case, Theorem 2.2 can be written as follows: 

Let b(u, v) = (1 — u)a(v) — u/3(v), and (Vn,Un) solution of the stochastic 

model (Sn)- Then, the process R N {t) = VN (u N {t) - u N (0) - J* b(u N (s), V N (s))ds 

converges in law, as N — > oo, towards the process R(t) defined as a stochastic 
integral: 

R(t) = I y/(l ~ u(s))a(V(s)) + u(s)/3(V{s))dB s 
Jo 

where B is a standard brownian motion and u(t), V(t) is the unique solution of: 

V = f(V,u) 

u=(l- u)a{V) - u/3(V) 

VN, u(0) = u = u N (0) 

VA^, V(0) =V Q = V n (0) 

This result provides the following degenerate diffusion approximation (Vat, Un), 
for N sufficiently large: 

dV N (t) = f(V N (t),u N (t))dt 
du N (t) = [(1 - u N (t))a(V N (t)) - u N (t)f3(V N (t))~\ dt + a N (u N (t), V N {t))dB t 

o-at(w) 2 = [(1 - u)a{v) + uf3(v)} = ^A(u, u) 

Let g(u, v) — \n(u, v)[^^i + (u, v) + -^^(u, v)] = (1 — u)a(v) + uf3(v). 
Note that in the multidimensional case, the real valuedfunction g above becomes 
a d x d-matrix. Since the different channel types j are supposed to independent, 
this matrix would be bloc diagonal, with blocs of size rj, thus assuring the 
independence of the q (rj-dimensional Brownian motions) in Theorem 2. 
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The blocs of size Tj are given by the matrix G^' of theorem 2, and arise from 
the calculation of the covariances: 



Gyj(x) = N\ N (x) J ZiZj/J,N{x,z) 

Proof of Theorem 2.2 We adapt the proof given by Kurtz [25]: we prove the 
convergence of characteristic functions plus tightness. The tightness property 
follows from the inequality: 

tN 

P[sup \R N (s)\ > 5} <-^|| 5 |U 

Let <j) N {t,9) = E[e ieR *V>] the characteristic function of R N . Let h(M N (t)) = 
e^W, VNM N (t) = R N (t), tP(u) = ^""^WV 2 , =e iu_ x _ iu = 
u 2 ip(u) — u 2 /2. We then have: 

<f> N (t,0)-l = E[h(M N (t))} - h(0) 

E[Aat(s) f h(w - u N (s) + M N (s)) - h(M N (s)) 
(w — ujv(s))/i'(Mjv(s))/ijv(s, dto)]ds 

E[e iei?N(s) Ajv(s) / C(6 , ViV(w-u J v(s))/ijv(s,dw)]ds 



E 



l e **xl«)\ N ( a ) { N9 2 (w-u N {s)) 2 fi N (s : dw) 



ds 



E[e ieRN ^X N (s) / N9 2 (w~u N (s)) 2 



En 



x tp (t/N9(w — ujv(s))^ (j>n(s, dw)]ds 

where Xn{s) stands for Ajv(wat(s), Vn(s)) and /zjv(s, rfxu) for /zjv(mat(s), Vat(s), dw). 
The second term in the last equality, call it K~n(9), converges to as N —> oo by 

dominated convergence, and because ip (yN9(w - u N (s)fj = ip (±9/Vn) -> 



as hmi/'fu) = 0. So we have: 

u— fO 



Mt,0)-i 



i 

2 7o 



2 6 



i0R N (s)a2 



g(u N (s),V N (s)) 



ds + K N {9) 



9 2 g(u{s),V( S ))0 N (s,9)ds 



1 



2 Jo 



6> 2 E 



Again, the second term in the last equality, call it Jn(9), converge to as 
N — >• oo, because of the convergence of ujv and Vjv to u and y.(cf. Theorem 



32 



2.1) 

By Gronwall lemma, we conclude that 4>N(t 7 8) — > <j)(t,9) with: 
4>(t,6)=cxp{-±9 2 J g(u(s),V(s))ds} ® 

5.2 Functional central limit theorem 

Let (Vn,un) be the solution of the simplified stochastic model (SW) and (V, u) 
of the deterministic model (D) introduced in the Example of section 2. Consider 
the process: 

( P N \ ( VN(u N -u) \ 
\Y N J \VN(V N -V)J 

If the initial conditions satisfy (P N (0),Y N (0)) = (0,0), the 2-dimensional pro- 
cess (Pn,Yn) converges in law, as N — > oo, towards the process (P,Y), with 
characteristic function: 



E 



j(e 1 p(t)+0 2 Y(t)) 



= e eiA(t)+e 2 2 B{t)+9 1 e 2 c{t) 



The functions A, B and C are solutions of the system: 



A' 

B' \ = 

a 



2b' u b' v \ ( A 
V' v f'u 5 I + 



iA(V,«) 





(M) 



with initial conditions (0,0,0), and with X(v,u) — — u)a(v) + u/3(v). 
Proof of Theorem 2.3 Just as in the proof of Theorem 2.2, let us define: 



E 



i(6 1 P N (t)+9 2 Y N (t)) 



Let us also define Zn = {un — u,Vn — V), Xjy = (un, Vat), X = (u, V), and 
h(x,y) = e l ^ R{6lX+e2V \ Then: 

<f>(t,6)-l = E[h(Z N (t))-h{Z N (0))} 

= ( E[NX(X N (s)) ( {h(w-u(s),V N (s)-V(s)) 

JO J E N 

- h(Z N (s))}fj,(X N (s),dw) 

- b{X{s))h' x (Z N (s)) + (/ (X N (s)) - f (X(s))) h' y {Z N { S ))]ds 
So<t> N (t,9)-l = G N (6,t) + H N {6, t) with 



G N (9,t) = J v[n N (s){(e ie ^ N - 1)^+ + ( e -^v0v/jv _i) M _}] 
n N (s) = NX(X N (s))h(Z N (s)), =n +/ _(X N (s)) 



H N (9,t) 



[ m 

Jo 



Nb{X{s))h{Z N {s)) 



ds 



N {/ (X N (s)) - f (X(s))} 
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Then in order to use the asymptotic development of e x when i^Owc introduce 
the function K{u) = e m — 1 — iu + u 2 /2. Then, knowing that /U + + fi- = 1: 

G N (6,t) = J* E Q N (a) {*^=(m+ - /OOM*)) - H + A-((9i/>/F)|l d.s 
Since fe(x) = A(x)(/i + (x) — we have : 



GW(M) 



/* 



Therefore: 

0jv(t,0)-l 



+ 
+ 

E 



E 



E 



i6WNb(X N (s))h(Z N (s)) 



ds 



~±\(X N (8))h(Z N (8)) 



ds 



[ E,\NK(9 1 /y/N)h(Z N (s) 
Jo L 



--6 2 2 \{X N {s))h{Z N {s)) 



ds 



ds (A) 

hiZffWWVN {b(X N (s)) - b(X(s))}] ds (B) 

+ f E \h(Z N (s))i9 2 VN{f(X N (s)) - f(X(s))}] ds (C) 
Jo 1 J 



+ 



E 



+ 



E 



h{Z N {s))NK{ei/yjN)\{X N {s)) ds (£>) 



Using the derivatives of b and /, and the convergence of Xjy to X we can make 
a development of the sum B + C: 

B + C = f E[h(Z N )iVN{(u N -u)(6 1 b' u + 6 2 fL) 
Jo 

+ (V N - + 9 2 f' v )}]ds} + e N (t, 9) 

where we dropped the s and where b' u , b' v , f' u , f' v are taken at Ajv(s). 

Noting that h{Z N )i\/N{u N -u) = h' x (Z N ) and h(Z N )i\/N(V N -V) = h' y {Z N ), 

we have: 

B + C = f E [h' x {Z N )(6tf» + 9 2 f u ) + ti v {Z N ){etf v + 9 2 .Q] 
Jo 

And the term D converges to zero as N — > oo by dominated convergence since 
K(u)/u 2 is bounded and converges to 0. 

As we have the convergence in Theorem 2.1 of A at to A, we get the convergence 
of <pN(t,0) to ^(f, 6), satisfying: 



1 



(t,6) = --9l\(X{t))^{t,9) + {9 1 b' u {X{t))+9 2 f' u {X{t))) 



(6 1 b' v (X(t)) + 6 2 f>(X(t))) 



Wo 



M 



Tightness stems from the Markov property and the following estimate ob- 
tained in the proof Theorem 2.1: 



sap J\VN[V N (t) - v(t)] \\ 2 + J2 

3 = 1 



a<t<T 



3 0) 



< exp 



(*) 



(*) 



(A/N)Ne 



> A 



— B(T)T 



CT 



The announced convergence in law follows. 

To solve the PDE, we set tf(i,0) = e 9 ? A ( t )+ fl i^ c ( t )+ 9 2 B W. Then, substitut- 
ing in the initial equation, and identifying the coefficients, we get the system 
(M). ® 

Proof of Theorem 2.4 We want to prove that the process Z has the same 
law as the limit as N — > oo of the difference between the Langevin approxima- 
tion linearized around the deterministic solutions and the deterministic solution 
itself, scaled by >/N. We write it in the general case, not only in dimension two 
as above. First we identify the equations satisfy by the moments of Z starting 
from the equation satisfied by the characteristic function. We make the ersatz: 

The matrix T t corresponds to the variance/covariance matrix. We plug this 
expression into the equation satisfied by ip as given in theorem 2.3: 



k dxi d6i 2 ^ k 
j=i {leL fe=i ' 1 fe,i=i 



+ EE* 

m=l leL 



dxi dOi 



The ensemble of indices L can be writen L = L v U L u where L v = {1 < m < p} 
and L u = {(j, k), 1 < j < q, 1 < k < rj}. To identify the equations satisfied 
by Tab we distinguish the following cases: 



a E L v and b E L v : \V ab = 



leL 



Oil] 
dxi 



df 



a e L„ and 6 e L„, & = (j, k): \T' ab = ^ 

leL 



dx L 
dxi 



J- a/ + 



a e L„, a = (j, fc) and & E L v : \T' ab 



E 



Ti — 1 6 

ox. 
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a G L u , a = (j, k) and b € L Vl b — (j' , k'): 

dbj^k , <96 7 < jfc < 



-r' =V 

2 oo / y 



r&z h — ?h — r a( 



2 k,k' L 3=J 



We then write the equations satisfied by K^ N \t) = VN(X N {t) - X(t)) = 
(Yff, Pjj ), where X N is the Langevin approximation defined in section 2.4, 
and where X is the deterministic limit: 



dYS 1 = VN(f m (X N ) - f m (X))dt 

dP^ k = VN{b j>k {X N ) - b j>k (X))dt + *W{X N )dWi 

When we linearize around the deterministic solution, we obtain the following 
equations: 



/'"„-<»>, 



leL 1 



where the terms -^Lfij^, comes from the linearization of aj? k ,(X N ), we do 
not need to specify them here because they go to zero as N — > 00. 

It is now clear that the moments equations for this linear diffusion system 
tends the system satisfied by T a b as N — > 00. 

Proof of Theorem 2.5 The convergence of Xm to X a.s. uniformly on finite 
time intervals, obtained in Theorem 2.1, implies that tat — > r a.s. In order apply 
Theorem 2.3, let us introduce Zn through the following decomposition: 

X(t n ) + ^=Z n (t n )^ -<t>{X{r N )) 



N{<I>(X(t))-<(>{X{t n ))) = Vn 

- VN<f>(X N (r N )) 



As N — > 00, we claim that the right hand side converges in law to V^>(X(r)).Z(r) 
since VN <j)(X n (tn)) converges in law to zero. Indeed, as <P(Xn{tn)) < and 
<I>(X n (t n )) > 0, 



'N<j>(X N ( TN ))\ < \VN(cP{X N (T N ))-(j,{X N {T^))\ 

There exists 8m on the line between Xn(tn) and Xn{t^) such that 

\Vn(4>(X n (t n )) - </>{X n {t„))\ = \V<j>(6 N ).(Z N (r N ) - Z n {t n ))\ 

which converges in law to zero since Zn — > Z and Z is continuous. The claim 
follows. By continuity, </>{X(t)) = 0, so that \/N(<j>(X (t)) - <j>(X(T N ))) is 
asymptotic to 

-V^(X(t)).F(X(t))Vn(t n -t) 
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Thus y/N (rjv — t) converges in law to 7r(r). To finish the proof we remark that 
\/N(X n (t n ) -X(t)) = Z n {t n ) + VN(X(t n ) -X{t)) which converges in law 
to Z(t)+tt(t)F(X(t)). 

A Comparison between two deterministic limits 
of different stochastic Hodgkin-Huxley mod- 
els 

We want to compare the two following systems deterministic (A.l) and (A. 2), 
with /, a,/3 continuously diffcrcntiable functions, a and j3 non-negative, k an 
integer > 1: 



~-f{V,u k ) 
w - (1 - u)a{V) - up(V) 



dt (12) 



du 



' f = f<y, Xk ) 

*l± = (k-j + l)x j _ 1 a(V) + (j + l)x j+1 l3(V) 
-x 3 (jp{V) + {k-j)a{V)) 
k VO < j < k 

System (A.l) corresponds to the classical "Hodgkin-Huxley" model, with 
only two variables for simplicity, and the system (A. 2) is a (fc + 2)-dimensional 
system, where Xj, < j < k is the proportion of channels in the state j, and 
j = k is the open state. 

Proposition Let Vq £ R and u e [0, 1]. If the following conditions on the 
initial values are satisfied: 

V(0) = V(0) = V Q and VO < j < k, C j k u{Q) k - ] {l - u(0)) j = x k -j(0) = 
C J k Uo~ j (l - u o y Then, for all t > 0, V(t) = V(t) (same potential) and u(t) k = 
Xi~(t) (the proportion of open channels is u{t) k ). 

Moreover, for all 1 < j < k, for all t > 0, x k -j(t) = C J k u(t) k - j (l - u(t)) j 

Proof Consider (V, u) the unique solution of (1) for V(0) = Vq and u(0) = Mo- 
Let y 3 (t) = C J k u(t) k -i(l - u(t)) j , < j < k. Then {V,y k , ...y Q ) is a solution of 
(2) (just need to compute y'j and write it in function of yj-i and yj+i)- As the 
initial values are equal (by hypothesis) : x k ^j(0) — C^u 1 ^ \1 — u n y = y k ^j(0), 
by uniqueness (V, y k , j/o) = {V, x k , ...xo) for all t > 0. 

Remark The result is essentially the same for more complicated Markov 
schemes, as the sodium multistate Markov model. 
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B Moments equations for linearized Langevin 
approximation 

From Theorem 2.2, one can build a diffusion approximation (Vjv,ejv) of the 
stochastic hybrid process (V/v, ejv) given in the Example of section 2.1: 



f dVi v (t) = /(Vi v (t),e JV (t))dt 



de^t) = b(e N (t), V N {t))dt + V ^"f'"'' ^ 
6(u,e) = [(1 - e)a(w) - e/3(»] 
[ A(v,e) = [(1 - e)a(v) + ep(v)} 

We want to write the moments equations for the linearized version of 



Pn 
Y N 



VN {e N - e) 

Vn (y N - v 



with (V 1 e) the deterministic solution. The linearized equations are given by: 
dYk = (/{rltf + f' e Pk)dt 



{ dP^ = {b' v Y^ + b' e Pk)dt 



2^NX t i^V Y N + K P k) 



dB t 



with A t = X(V(t),e(t)). We define mf 
mO 2 ], = E[(i$ - m 2 ) 2 ] and C£ 
the following system of 5 equations: 



■ JVJ) 



E[P£], Sf 



AT 



E[(yf 



AT 



mi ) 2 ], 5^ = E[(P£ - m 2 ) 2 ] and = E[(Y^ - mi)(P£ - m 2 )]. Then we have 



dmf 

-f dt » 

L dm^ 
dt 



6' y mi + 6' e m 2 



dSf 

dt 
dSf 

dt 



dC£ 
dt 



(A'^mi + A' e m 2 ^ 



2/{,Si + 2/^C 12 
26^2 + 26^12+ [VA7 + 

(i) 2 *+(A) 2 ^+ 2 &- 

6V5 1 + /^ 2 + (/(, + 6' e )C 12 



At the limit N — > oo and with A = — 26*2, P = — 2S 1 ! and C = — Ci 2 this system 
is the same as the one found in application of Theorem 2.3 in section 3. 
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C Moments equations for the Morris-Lecar sys- 
tem 



The moments equations used in section 3.5.1 and 3.5.2 are the following linear 
non-homogeneous system of differential equations: 







$711 




Bi 


S n 




Sn 




B 2 


Sy 


= M(t) 


Sn 
Cmv 


+ 










C n v 



















with 



M(t) 



! c)F m 
dm 

o 
o 

dF v 
dm 








o dF„ 
Z dn 




dn 








dW m 



dV 




dFv 
dV 



o dF m 
Z dV 



2fi 



dm 



dF n 
dV 





odF n 
Z dV 
02 9Ji 



dFy 
dV 



dn 



dn 



+ 



dF„ 
dV 



dF m 
dm 







dFv 

Ml 

dm 

+ 



dFs. 

dn 



all the functions being evaluated at X(t) = (V (t) , m(t) , n{t)) solution of (3.5- 
3.7) and with B x {t) = (1 - m(t))a m (V(t)) + m(t)(3 m (V(t)), B 2 (t) = (1 - 
n(t))a n (V(t))+n(t)p n (V(t)). 
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